% Scripts to test the difference between day and night (or all) data for
% tidal harmonics estimation. The station chosen is ASC because it has both
% M2 and N2 signals


% Get a list of stations

S = dir('/data/backup/mnair/obs_mag_data/update_2011/secular_removed/*.nc');
for i = 1:length(S),
    robstn(i,:) = S(i).name(1:3);
end;
a = (datenum(2011,1,1)- datenum(1995,1,1))*24*60;
t = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2010,12,31,23,59,30));
periods = [4 4.8 6 8 11.967236 12  12.4206012 12.6583 23.934472 24 25.891];

%periods = 12.4206012;

%periods = 12.421;

x = [cos(2*pi*t/( periods(1)/24))'       sin(2*pi*t/( periods(1)/24))'];

for j = 2:length(periods),
    
    x = [x cos(2*pi*t/( periods(j)/24))'       sin(2*pi*t/( periods(j)/24))'];
    
end;


%nstn = 11; % ASC station


[dummy,nstn,ib] = intersect(robstn,'MBO','rows')
    
ncid = netcdf.open(['/data/backup/mnair/obs_mag_data/update_2011/secular_removed/' ... 
    S(nstn).name],'NOWRITE');


X_ID = netcdf.inqVarID(ncid,'Magnetic_Field_X');
Y_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Y');
Z_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Z');
z_data = netcdf.getVar(ncid,Z_ID,0,a);

latitude  = double(netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'geospatial_lat'));
longitude = double(netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'geospatial_lon'));

z_data = double(z_data)/10;
z_data(z_data == 99999.9 ) = NaN;
L = isnan(z_data);
k = find(L==0);
z_data(L) = 0;

netcdf.close(ncid);


%Found the the data for periods > 2005 (low geomagnetic activity)

LL = t > datenum(2005,01,01);

x(~LL,:) = [];
t(~LL) = [];
z_data(~LL) = [];
L(~LL) = [];


%robust fit
%all data

spectra_all = robustfit(x(~L,:), z_data(~L));

% calculate the amplitudes

spectra_data_all_rob = abs(complex(spectra_all(2:2:end),spectra_all(3:2:end)));


%least square solution

spectra_all_ls = x(~L,:)\z_data(~L);

spectra_data_all_ls = abs(complex(spectra_all_ls(1:2:end),spectra_all_ls(2:2:end)));


% only nightime data

%Get the LT
t = t + (longitude/15)/24; %shift the axis to the local time
L = L | (t-floor(t)  < 18/24 & t-floor(t)  > 6/24)';

%robust fit of the night data
spectra_night_rob = robustfit(x(~L,:), z_data(~L));
spectra_data_night_rob = abs(complex(spectra_night_rob(2:2:end),spectra_night_rob(3:2:end)));


%ls fit of the night data
spectra_night_ls = x(~L,:)\z_data(~L);
spectra_data_night_ls = abs(complex(spectra_night_ls(1:2:end),spectra_night_ls(2:2:end)));

% spectra using only M2 component
x = [cos(2*pi*t/( periods(7)/24))'       sin(2*pi*t/( periods(7)/24))'];
%ls fit of the night data for M2 component
spectra_night_ls_m2 = x(~L,:)\z_data(~L);
spectra_data_night_ls_m2 = abs(complex(spectra_night_ls_m2(1:2:end),spectra_night_ls_m2(2:2:end)));

% sweeping in frequency for M2 component LS
prd = 12.415;
clear spectra_data_ls_m2_n
for i = 1:100,
    
x_m2 = [ones([1,length(t)])' cos(2*pi*t/( prd/24))'  sin(2*pi*t/(prd/24))'    ... %M2
    ];  
spectra_night_ls_m2 = x_m2(~L,:)\z_data(~L);
spectra_data_ls_m2_n(i) = abs(complex(spectra_night_ls_m2(2),spectra_night_ls_m2(3)));
prd = prd + 0.0001;
end;




% sweeping in frequency for M2 component ROB TOO SLOW
prd = 12.415;
clear spectra_data_rob_m2_n

for i = 1:100,
    
x_m2 = [ cos(2*pi*t/( prd/24))'  sin(2*pi*t/(prd/24))'    ... %M2
    ];  

spectra_night_rob_m2 = robustfit(x_m2(~L,:),z_data(~L));
spectra_data_rob_m2_n(i) = abs(complex(spectra_night_rob_m2(2),spectra_night_rob_m2(3)));
prd = prd + 0.0001;
end;

